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Abstract 

In this work, the sharp interface limit of the degenerate Cahn-Hilliard equation 
(in two space dimensions) with a polynomial double well free energy and a quadratic 
mobility is derived via a matched asymptotic analysis involving exponentially large 
and small terms and multiple inner layers. In contrast to some results found in the 
literature, our analysis reveals that the interface motion is driven by a combination of 
surface diffusion flux proportional to the surface Laplacian of the interface curvature 
and an additional contribution from nonlinear, porous-medium type bulk diffusion, 
For higher degenerate mobilities, bulk diffusion is subdominant. The sharp interface 
models are corroborated by comparing relaxation rates of perturbations to a radially 
symmetric stationary state with those obtained by the phase field model. 


1 Introduction 

Phase field models are a common framework to describe the mesoscale kinetics of phase 
separation and pattern-forming processes ga El]. Since phase field models replace a 
sharp interface by a diffuse order parameter profile, they avoid numerical interface track¬ 
ing, and are versatile enough to capture topological changes. Although such models can 
be constructed starting from a systematic coarse-graining of the microscopic Hamilto¬ 
nian | dll 30, (29J {28], the use as a numerical tool to approximate a specific free boundary 
problem requires in the first instance careful consideration of their asymptotic long-time 
sharp interface limits. 

In this paper, we will mainly focus on the Cahn-Hilliard equation for a single con¬ 
served order parameter u = u(x, t), 

«t = -V-j, j = -M(u)V/r /i = —£ 2 V 2 u + (la) 

with a double well potential 

f{u) = (1 -u 2 ) 2 /2 (lb) 
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and the degenerate, quadratic mobility 

M(u) = (l-u 2 ) + , (lc) 

on a bounded two-dimensional domain with boundary conditions 

Vu • n = 0, j • n = 0 (Id) 

at Here, (•)+ is the positive part of the quantity in the brackets, x represents the 
two-dimensional spatial coordinates, t is the time, \i the chemical potential, j the flux, 
and n the outward pointing normal to <917. Boldface characters generally represent two- 
dimensional vectors. Both the potential and the mobility are defined for all u. The 
mobility is continuous but not differentiable at u = ±1. 

The case of a Cahn-Hilliard equation with a constant mobility has been intensively 
discussed in the literature. In particular, the sharp interface limit e —> 0 was determined 
by Pego [46] . and subsequently proven rigorously by Alikakos et al. [3j. On a long time 
scale, t = the result is the Mullins-Sekerka problem jll| . In particular, the 

motion of the interface between the two phases is driven by flux from bulk diffusion. 

In contrast, Cahn-Hilliard equations with degenerate mobility are commonly expected 
to approximate interface motion by surface diffusion [43] on the time scale t = 0(e -2 ), 
where the interface velocity v n is proportional to the surface Laplacian A s of the interface 
curvature k, 

v n oc A s n. (2) 

We note that the surface Laplacian is equal to d ss K in two space dimensions, where s is 
the arclength. In fact, for the case of the degenerate mobility M(u ) = 1 — u 2 and either 
the logarithmic free energy 

/(«) = \ e K 1 + u) ln(l + u) + (1 - u) ln(l - it)] + 1 - u 2 ), 

with temperature 9 = 0(e a ), or the double obstacle potential 

f(u) = 1 — u 2 for |u| < 1, f{u) = oo otherwise, 

Cahn et al. [ T8] showed via asymptotic expansions that the sharp interface limit is indeed 
interface motion by surface diffusion (f2j). 

Although the logarithmic potential and the double obstacle potential as its deep 
quench limit are well motivated, in particular for binary alloys, ngna Sana ehi GEi Ha 
da, other combinations of potentials and mobility have been used in the literature as a 
basis for numerical approaches to surface diffusion [20]. Those models are often employed 
in more complex situations with additional physical effects, such as the electromigration 
in metals [42], heteroepitaxial growth [ 49] - anisotropic fields mm, phase separation 
of polymer mixtures [5811571 and more recently in solid-solid dewetting [ 34] and coupled 
fluid flows El ED El- In those models, a smooth polynomial double-well free energy is used 
in combination with the mobility M (u) = 1 — u 2 or the degenerate biquadratic mobility 
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M(u) = (1 — it 2 ) 2 for |u| < 1. A smooth free energy is numerically more convenient to 
implement, especially in a multiphyscial model, as it avoids the singularity present in 
either the logarithmic or double obstacle potential. Authors typically justify their choice 
of mobility and free energy by adapting the asymptotic analysis by Pego [46] and Cahn 
et al. |18j to obtain the interface motion ([2]) in the sharp interface limit. 

Interestingly, Gugenberger et al. [33]. recently revisited some of these models and 
pointed out an apparent inconsistency that appears in the asymptotic derivations except 
when the interface is flat. Other evidence suggests that the inconsistency may not be 
a mere technicality but that some bulk diffusion is present and enters the interfacial 
mass flux at the same order as surface diffusion. This was observed for example by 
Bray and Emmott [T5] when considering the coarsening rates for dilute mixtures, and by 
Dai and Du [22] where the mobility is degenerate on one but is constant on the other 
side of the interface; the papers by Glasner [32] and Lu et al. m also use a one-sided 
degenerate mobility but consider a time regime where all contributions from the side 
with the degeneracy are dominated by bulk diffusion from the other.) In fact, an early 
publication by Cahn and Taylor m remarked that using a biquadratic potential might 
not drive the order parameter close enough towards ±1 to sufficiently suppress bulk 
diffusion, citing unpublished numerical results. Diffuse interface models for binary fluids 
with a double well potential and a quadratic mobility M(u ) = 1 — u 2 or M(u) = (1 — ri 2 )+ 
are investigated in mm- However, in both studies, the leading order expressions for 
the interface motion do not contain bulk diffusion contributions. 

In this paper, we aim to resolve the apparent conundrum in the literature, and re¬ 
visit the sharp interface limit for ©■ We will obtain a sharp interface model where 
the interface motion is driven by surface diffusion, i.e. the surface Laplacian, and a flux 
contribution due to nonlinear bulk diffusion either from one or both sides of the inter¬ 
face, depending on the nature of the solutions for u in the outer regime. The matched 
asymptotic analysis is rather subtle, and involves the matching of exponentially large 
and small terms and multiple inner layers. 

The paper is organised as follows: Section 2 approximates solutions of m which 
satisfy |u| < 1; Section 3 considers the asymptotic structure of the radially symmetric 
stationary state, which demonstrates the matched asymptotic expansion and exponential 
matching technique in a simpler setting; Section 4 returns to the general 2D time depen¬ 
dent problem; Section 5 briefly discusses the sharp interface limit for a class of solutions 
with the mobility M(u ) = |1 — v?\ where |u| < 1 is not satisfied, and for the Cahn- 
Hilliard model with a biquadratic degenerate mobility M(u) = ((1 — rt 2 )+) 2 ; Section 6 
summarises and concludes the work. 

2 Preliminaries 

In this paper, we are interested in the behaviour of solutions to (lla[) describing a system 
that has separated into regions where u is close to ±1, except for inner layers of width e 
between them, and evolve on the typical time for surface diffusion, t = 0(e~ 2 ). We thus 
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(3a) 


rescale time via t = e 2 t, so that the Cahn-Hilliard equation reads 

e 2 <9 T ii = V-j, j = M(u)V/x, fj, = — e 2 V 2 u + f'(u), 

and we keep the boundary conditions on 517, 

Vu • n = 0, j • n = 0. at 50. (3b) 

We will denote the subsets where u > 0 and u < 0 by 0 + and 0_, respectively, and 
identify the location of the interface with u = 0. Moreover, we assume that 0 + is convex 
unless otherwise stated, and has 0(1) curvature everywhere. We will focus on solutions 
of ([3k, b) that satisfy |tt| < 1. The existence of such solutions has been shown by Elliott 
and Garcke [23] . 

The general procedure to obtain a description of the interface evolution is then to 
consider and match expansions of (J3^i,b) , the so-called outer expansions, with inner ex¬ 
pansions using appropriate scaled coordinates local to the interface. The approach as¬ 
sumes that the solution of ([3k, b) is quasi-stationary i.e. close to an equilibrium state. 
Unfortunately, it is not obvious what the appropriate nearby equilibrium state could be 
in the situation we consider here. The problem arises because equilibrium solution to 
([3^, b) with constant y does not generally satisfy the bound |u| < 1 inside of f2_|_ [46J. 

It is helpful to revisit the standard matched asymptotics procedure for ([3k, b) to 
understand the implications of this observation. Notice that the time derivatives drop 
out of the lower order outer and inner problems. The leading order inner solution for 
the double well potential is simply a tanh-profile, which matches with ±1 in the outer 
solution; the corresponding leading order chemical potential is zero. To next order, 
the inner chemical potential is proportional to n, and this supplies boundary conditions 
for the chemical potential in the outer problem via matching to be y,\ = c\K. Here, //i 
denotes the first non-trivial contribution to the chemical potential in the outer expansion, 
fjL = Efi\ +0(£ 2 ), and ci represents a fixed numerical value. It is obtained from a detailed 
calculation along the lines of section [31 which in fact shows that ci > 0. It is easy to see 
from the third equation in (I3a[) that the outer correction u\ for u = ±1 + eu\ is given 
by u\ = fii /thus u = ±1 + cine/A + 0(e 2 ) near the interface. Inside f2 + , we 
therefore have that the outer solution u > 1. Notice that we have used that / is smooth 
at u = ±1 — for the double obstacle potential, there is no correction to u = ±1 in the 
outer problem, see |18j . 

The resolution to the above conundrum comes from the observation that for a degen¬ 
erate mobility, slowly evolving solutions can arise from situations other than constant 
H once |rt| gets close to 1. To obtain an indication of how such solutions evolve, we 
look at numerical solutions of the radially symmetric version of ([3k, b) on the domain 
17 = r < 1}, where r = (x 2 + y 2 ) 1//2 , starting with a tanh as initial profile such 

that Uinit(r’) < 1. The spectral method we used is briefly described in the appendix. The 
numerical solution at a later stage as shown in Fig. [T| is positive for r < 0.5 and negative 
for r > 0.5. Notice while for r > 0.6 the solution for u levels out into a flat state that is 
larger than —1 by an amount of O(e), for r < 0.4 the solution is much closer to u = 1. 
Closer inspection shows that u has a maximum which approaches u = 1, say at r = r*. 
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Figure 1: The long-time solution u for the radially symmetric degenerate Cahn-Hilliard 
equation m for different initial data and different mobilities. In (a, left panel), the 
mobility is (TTcl) and initial data is bounded within [—1,1], while in (b, right panel) it 
exceeds 1 and —1 to the left and right, and the mobility is replaced by M(u ) = |1 — it 2 1, 
respectively. In both panels, the initial data is shown by dashed lines while the long-time 
solutions for e = 0.05 are given by solid lines and have converged close to a stationary 
state. In (a), this stationary profile is bounded between [—1,1], where we emphasize that 
u in the left inset is still below 1 (dashed line in the inset), while in (b), the upper bound 
1 is exceeded for r less than about 0.4 (see left inset in (b)). Notice that in both (a) and 
(b), the value for u for r > 0.7 is close to but visibly larger than — 1, by an amount that 
is consistent with the 0(e) correction predicted by the asymptotic analysis (for (a) in 

csD). 


The maximum of u may touch u = 1 in either finite or infinite time. In either case, the 
solution in splits into two parts to the left and right of r*. The flux between the 
two parts is very small, and this suggests that they are nearly isolated from each other. 
In particular, they do not have to be at the same chemical potential. Since we are only 
interested in the phase field where it determines the evolution of the interface, we cut 
off the part with r < r*, and consider the remaining part r > r* as a free boundary 
problem. 

Returning to the general case of not necessarily radially symmetric solutions, we 
introduce a free boundary T near the interface inside Q+, and cut off the parts of the 
solution further inside of f2 + . At T, we impose 

u = 1, np • j = 0, np • Vit = 0. (3c) 

Notice that in addition to u = 1 and vanishing normal flux, a third condition has been 
introduced at T. This is expected for non-degenerate fourth order problems and permits 
a local expansion satisfying ll3cl) that has the required number of two degrees of free¬ 
dom [35] . Indeed, expanding the solution to 0 in a travelling wave frame local to T 
with respect to the coordinate r? normal to T gives u = 1 — arf + 0(r/ 3 ), where a and the 
position of the free boundary implicit in the travelling wave transformation represent the 
two degrees of freedom. 

Also observe that if it > — 1 by 0(e) as suggested by the numerical solution in 
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Fig. Oja), then M(u ) = 0(e). Since y = 0(s), we expect a nonlinear bulk flux of order 
0(e 2 ) at the interface arising from I2_. This is the same order as the expected flux from 
surface diffusion. Indeed, as shown below, both contributions are present in the leading 
order sharp interface model (158d[) . 

Another scenario is conceivable if the mobility is changed to 11 — v 2 1. Then, with an 
appropriate initial condition, we obtained numerical results for the radially symmetric 
case which suggest a solution that is not confined to |n( < 1 and which in fact converges 
to the usual stationary Cahn-Hilliard solution (considered, for example, in [45 j ) for which 
y is constant in 17, and u is larger than one in most of These results are shown in 
Fig. IHb). In this case, bulk fluxes from both and contribute to the leading order 
interface dynamics, see section l5Tl 


3 Radially symmetric stationary solution 

By setting u T = 0 in Q for a radially symmetric domain 0 = {(x,y);r < 1} and radially 
symmetric u = u(r), where r = (x 2 + y 2 ) 1 ' 2 and then integrating we obtain 

£ 2 d / du\ , 9 . ., . 

-— ( r— + j? - 2 u(u 2 - 1) = 0, (4a) 

r dr \ dr J 

u'{ 1) = 0, (4b) 

u(r*) = 1, u'{r*) = 0. (4c) 

The point r* represents the location of the free boundary T that needs to be determined 

as part of the problem. The chemical potential y is constant that needs to be determined 

by fixing the size of the 0_j_. This can be done by specifying the or, simpler, the 

position ro of the interface, 

u(r 0 ) = 0. (4d) 

Note that if we do not consider a free boundary T and impose u'( 0) = 0 instead of (I4cl) . 
then there exist exactly two solutions (which can be discerned by the sign of rt( 0 )) as 
was shown in [45] ■ We will now investigate (j4[) in the sharp interface limit e —> 0 using 
matched asymptotics. There is one outer region away from the interface, and two inner 
layers, one located at the interface ro and one located at r*. 

Outer region 

Inserting the ansatz 


u = u 0 + £U\ H-, V = Vo + erji-\ -, 

into Mali and (]4b]) and taking into account that the chemical potential 77 is a constant 
quickly reveals that uo, u\ and 112 are also constants. Their values are fixed by standard 
matching, that is, they are equal to the limits of the inner solutions as p —> 00 , which 
therefore have to be bounded in this limit. 
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Inner layer about the interface 


To elucidate the asymptotic structure of the interface, we strain the coordinates about 
ro and write 

r ~ r 0 ,-x 

P= ——, (5) 


so that for U(p) = u(r), and with the interface curvature k = 1/ro, we have 

U' 


U" + £ 


k 1 + ep 


+ T} - 2(U 6 - U) = 0, 17(0) = 0. 


Expanding U = Uq + eU\ + • • • , we have, to leading order, 

US ~ 2(t / 0 3 - Uq) = vo, Uq{ 0) = 0. 


( 6 ) 


(7) 


To match with the outer and the solution near T, Uq needs to be bounded for p —> Too, 
which gives 

Uq = — tanh p, vo = 0. 


To 0(e) we have 

U'{ - 2(3Uq - l)t/i = -m - kWq, C71(0) = 0, 
for which the solution that is bounded as p —> oo is given by 


( 8 ) 

(9) 


1 


1 


Ui = “ ——(771 + 2K)sech^p + — (3771 — 2At)sech z p — + - sinh 2p + — sinh 4p 


3 p 1 . 


1 


16 

H—(2 k — 771 ) 4- (2k — 3r?i)(2cosh2 p — 5 sech 2 p). 

8 48 


32 


( 10 ) 


Inner layer about T 

We centre the coordinates about the free boundary r = r* and write 

z = p + a, a = (r 0 r*) / e. (11) 

Substituting in the ansatz U = 1 + eU\ + e 2 U 2 + ..., we obtain, to 0(e), the problem 


with the solution 


U" ~ 4t7i = —pi, 

Ui(0) = 0, u[(o) = 0, 

Ui= 1 j(l- cosh22:) . 


(12a) 

(12b) 

(13) 
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Matching 

We first observe from (l4cl) that the location of the free boundary T in the inner coordinate 
p = —a satisfies U(—a) = 1, U'(—a) = 0. However, for e — > 0, we also have U —> Uq = 
— tanh(p) <1. To reconcile these conditions, we need to assume a —» oo as e —> 0. 
Matching of the inner expansions therefore involves exponential terms with large negative 
arguments p , or conversely for large positive z, which we deal with in the spirit of Langer 
|3TJ |. see also [38] . The solution centred at the interface is expanded at p —> — oo and 
the result written and re-expanded in terms of z = p + a. Notice that this change of 
variables can lead to terms changing their order in e if a has the appropriate magnitude. 
The solution for the layer around the free boundary T is directly expanded in terms of 
z oo and then the terms are matched between the two expansions. 

Expanding Uq and U\ for p —» — oo and substituting p = z — a gives 


U = 


1 - 2e -2fT e 2z 


+0(e 4z )) 


+ £ < 


^(2k — 3f7i)e 2<T e 2z + \(k - pi) 


B 



+0(s 2 ). 



D 


(z-a) 


e~ 2a e 2z 


+0(e 


4z\ 


(14) 


The inner expansion for U at z —> oo is 



—^-e~ 2z +0(s 2 ) 


G 


(15) 


Comparing terms in (I14p and (11511 of the same order of e functional dependence with re¬ 
spect to z, we notice first that the constant terms at 0(1) are already matched. Matching 
dO and [El yields 



(16) 


As a result, the term [B] is zero. Matching term [A] and [Fj we arrive at the condition 
2e _2fT = e/c/12, which we solve for a, giving 



(17) 


We can now determine the outer solutions. We note that in the more general, time 
dependent situation, the presence of a non-zero correction will give rise to a flux at 0(e 2 ). 
Using the limits of Uq and U± as p —> oo, we obtain 


u 0 = - 1 , 


K 


Ul= 6- 


(18) 
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Higher corrections 


At this stage, it is obvious that the matching is not yet complete to O(e), as the terms 
in (fT5l) and (I14|) . respectively, TTTI and iGl are non-zero and lack counterparts in the other 
expansion. This can be resolved by considering the next higher order solutions U 2 and 
U 2 , which, in fact, will also be useful in section [4j We include e 2 p 2 in the expansion for 
7 /, and allow for corrections to a via the expansion 


1, /24\ 

G = o lo S — + ecr i H-• 

2 V £K / 

The 0(e 2 ) problem at the interface is given by 

UZ - 2(3C/ 0 2 - l)U 2 = -m - KU[ + pk 2 Uq + 6U 0 Uf 

K? K 2 

= —P 2 -tanh° p — pK 2 sech 2 p-tanh p sech 2 p, 

6 3 

together with U 2 {0i) = 0 and boundedness for C /2 as p —> 00 . The solution is 

r /2 pn 2 1, ( 2n\ 1 9 / 23 0 9 9 

U-i = ----cosh2 p ( 772 - 1 - -pK 2 J + — sech 2 p f 5r/ 2 + ~ 2 P K 


(19) 


( 20 ) 


+ 


1 


pn 2 log ^e p ^ sech 2 p + ^-sech 2 p Li 2 (—e 2p ) 


-sinh 2 p (1 — 24 log cosh p) 

288 

—— tanhp (1 — 24 log cosh p -sech 2 p ) 4- ( — n 2 — 772 J sech 2 p 

96 \ 3 / 16 \ 6 J 

+ (^— (1 + 24 log 2) + r/2^1 sech 2 p (^ \ sinh 2p + sinh 4 p 

\3o J \ 8 4 32 

where Li 2 (x) is the dilogarithm function. 

For U 2 {z) we have 

c7" - 4C/ 2 + kC7( - 6C/ 2 + m = 0, 
t^2(0) = 0, C7'(0) = 0, 


( 21 ) 


(22a) 

(22b) 


which has the solution 

U 2 = (^) (cosh4z + 3e _2z (l + 4z) - 9) + e 2z 

+ (^j e_2z + ^(1 - cosh2z). (23) 

Expanding U = Uq + eU\ + e 2 U 2 + ■ • • for p —>• — 00 , substituting in p = z - cr and 
using csi leads to 


C/ 


1 - 


-( 

12 


,2z 


(1 — 2e<7i) + 


1 



+ £ 





(1 + 2e<Ji)e 2z + 



+ 0(e 3 ). 


(24) 
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Figure 2: Comparing the asymptotic and numerical results for (left) the position of the 
free boundary and (right) the chemical potential, for a range of e and r$ = 1/2. 


Similarly, the expansion for U = Uq + eU\ + e 2 U 2 + ■ ■ ■ as z —> oo is 


U = 1 + £— (1 — cosh 2 z) 

6 


+ £ 

+ 


1 (—Y~ 4z 

2 \ 12 / 


' e '" + 2 (ll) e- + (^)-(3e--(l + 4z)-9) 


12 


^2z 


+ (0 e 2z + ^(1 - cosh2z) 


(25) 


Now, we can match the e 2z at 0(e) and the e 2z at 0(e 2 ) terms, and arrive at, respec¬ 
tively, 


V ' 2 36 ’ 


3k 

ai = 16' 


(26) 


For completeness we note that the next order outer correction u 2 is again a constant 
equal to the limit of U 2 as p —> 00 , with the value 112 = 7k 2 /144. 

Figure [2] shows that the asymptotic results agree well with the position of T and 
the chemical potential obtained from numerical solutions of the ODE free boundary 
problem (H|), confirming the validity of the matched asymptotic results. The solutions 
were obtained by a shooting method with fixed r] using the Matlab package odel5s, with 
'u(l) and (l4cj) as the shooting parameter and condition. The value of rj is adjusted in an 
outer loop via the bisection method until ro = 1/2 is achieved to a 10 -10 accuracy. 


4 Sharp Interface Dynamics 

4.1 Outer variables 

Motivated by the stationary state, we now consider the asymptotic structure of the 
dynamical problem that arises for non-radially symmetric interface geometries. For the 
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outer expansions, we will use 

u = u 0 + eui + £ 2 u 2 H -, p = p 0 + + £ 2 H2 H -, j = jo +eji +£ 2 h H -■ 


4.2 Inner variables 

As in [46 . 33], we define the local coordinates relative to the position of the interface 
(parametrised by s), and write 


r(s,r,r) = R( 

where R, the position of the interface (, is 

u(R, 


s,t) + rn(s,r), 

(27) 

defined by 


t) = o, 

(28) 


and t = dR/ds is the unit tangent vector, and n is the unit outward normal. From the 
Serret-Frenet formulae in 2D we have that nt = <9n/<9s, thus 

|)l = n 0)> ^ = (l + r-«)t(s), (29) 

where t(s) is the unit tangent vector to the interface, and k is the curvature. We adopt 
the convention that the curvature is positively defined if the osculating circle lies inside 
The gradient operator in these curvilinear coordinates reads 


V = nd r + 


1 


-t d. 


1 + rn 

and the divergence operator of a vector field A = A r n + A s t reads 


V • A = 


1 


1 + rn 


d r ^(1 + + d s ^ ^ ^ rK^ s 


(30) 


(31) 


We let s and p = r/s be the inner coordinates at the interface, and let U(p,s,t), 
r/(p,s,T ) and J(p, s,r) denote the order parameter, chemical potential and flux written 
in these coordinates, respectively. In inner coordinates, the combination of the first two 
equations, in (|3a[) . and (|28p . become 

e 2 d r U - ev n dpU = V • (M(C/)Vr?), (32a) 

V = s 2 V 2 U + f \U), (32b) 

17(0) = 0, (32c) 
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with v n = R r • n. Using equations (I30[) and (1311) . we obtain 
V ■ (M(U)'V) = s~ 2 d p M(Uo)dp 

+£- 1 \d p {npM{U Q ) + M'{U Q )U-^dp - k P d p M(U 0 )dp\ 

+ | K 2 p 2 d p M(Uo)d p - npdp(^pM{U 0 ) + M\U Q )U^d p 
+d p ^pM , (Uo)U 1 + ^M"{U 0 )Ul + M’{U Q )U^d p 
+d s M(U 0 )d s \ +0(e). (32d) 


Notice that the corresponding expression for V 2 can be easily obtained from this by 
setting M = 1. 

Taking only the first equation in (I3a|) we have 


£ 2 d T U — £v n d p U = 


1 


1 + £pK 


£-'dp({l + £pK)J n )+d s ( Y ^J l 


(33) 


In inner coordinates, we will only need to know the normal component J n = n • J of the 
flux explicitly in terms of the order parameter and chemical potential. It is given by 

= £~ 1 M(Uo)dpT]o + M' (Uo)U\dpr}Q + M(U 0 )d p m 
+ £ (M(Uo)d p rj 2 + M'(Uo)Uid p r]i + M\U 0 )U 2 d pVo + ^M"(U 0 )Ufd p rio 


+ £ 


M(U 0 )d p r, 3 + M\U 0 )Uid p rj 2 + ( M'(U 0 )U 2 + -M"(U 0 )Uf d p r h 


+ M'(U 0 )U 3 + M {Uq)U\U 2 + -M" (Uq)U{ dp Vo 


+ 0(e 3 ), (34) 


which also motivates our ansatz for the expansion for J, given the obvious ansatz for the 
other variables, 


U = U 0 + £U\ + £ 2 u 2 H-, p = rjo + £7/1 + £ 2 r ]2 H-, 

J = £ ^J~i T Jq T £«Ji T e“J 2 + • • • . 


Moreover, we introduce 2 = p+a(s, t) as the coordinate for the inner layer about the 
the free boundary T, so that the order parameter, chemical potential and flux in these 
variables are given by U(z,s,r), fj(z,s,T ) and J(z, s,r), respectively, with expansions 

U = Uq + £U\ + £ 2 U 2 + • • • , fj = fjQ + £fj i + £ 2 ff 2 + • • • , 

J=£ 1 J_1 + Jq + eJl + J 2 + ■ ■ ■ ■ 
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Notice that the location where the two inner layers are centred depends on e and 
therefore, in principle, a and also R need to be expanded in terms of e as well. However, 
we are only interested in the leading order interface motion, so to keep the notation 
simple, we do not distinguish between a and R and their leading order contributions. 


We now solve and match the outer and inner problems order by order. 

4.3 Matching 
Leading order 

For the outer problem, we obtain to leading order 

V ■ jo = 0, j 0 = M(u 0 )V/u 0 , no = f'(u 0 ). (35) 

The requisite boundary conditions are V n Uo = 0, and n ■ jo = 0 on dQ. We have 

uq = -1, Ho = 0. (36) 

The leading order expansion about the interface reads, 

M(U 0 )d p p o = ai (s, r), f'(U 0 ) - d pp U 0 = Vo- (37) 

From the matching conditions, we require Uq to be bounded for p —> Too. In fact, 
U(p —> —oo) = —1, giving po 0- This implies ai = 0, therefore also po = 0, which we 
note matches with ^o- Moreover, from (I37 [) q and from (|34l) we have 

Uq = — tanhp, J n -i = 0. (38) 

The leading order approximation of the order parameter in the coordinates of the inner 
layer at T is easily found to be Uq = 1, and also for the chemical potential po = 0, and 


the normal component of the flux J n -1 = 0. 

O(e) correction 


The first two parts of the outer correction problem for (I3al) are automatically satisfied, 
since Ho = 0 and M(uq) = 0, by 



ji = 0. 

(39) 

The last part requires 

Hi = f"{uo)ui = Au\. 

(40) 

From (132 p. and noting that po 

= 0, we have 


d p (M(U 0 )d p p 1 ) = 0, 

Pi = -d PP Ui - Kd p U 0 + f"(U 0 )Ui, lh (0) = 0, 

(41) 


thus M(Uo)d p pi = J n fl is constant in p. Since J n> o has to match with jo, it is zero. 
Therefore, p\ = pi(s,t) does not depend on p. Now (1411) 9 and (141[) represent the same 
problem as ©. As such, the solution U±(p, s, r) that is bounded as p —> oo can be read 
off odd. 
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The 0(e) problem for the inner layer at T becomes 


m = -d zz U\+m x , 


( 42 ) 


with r/i that does not depend on z, supplemented with the conditions Ui(z, 0,r) = 1, 
U\ z (z, 0, r) = 0. This equation is the same as the 0(e) equation for the stationary state 
about the free boundary, and the solution is given by (USD- The inner layers about T and 
about the interface can be matched, as outlined in section [3j to obtain 



(43) 


0(e 2 ) correction 

Combining the first two equations in (IHaTl and expanding to 0(e 2 ) yields 

V • (M' (u 0 ) Ul Vin) = 0. (44) 

In view of the discontinuous derivative of M at u = uq = — 1, we remark that here and 
in the following we will use the convention that M'(± 1) denotes the one-sided limit for 
|u| —> 1 ~, in particular that M'(— 1) = 2, and likewise for higher derivatives. Equation 
(14011 provides a relation between p\ and u±. Thus, we have 

V • (/riV/ri) = 0 (45) 


with the boundary condition \7 n p\ 
(USD) at the interface, 


0 on dQ, and, from matching p\ with r /i (given in 
Mi = |k- (46) 


Expanding the second equation in (l3al) to 0(e 2 ) also gives us an expression for the normal 
flux 


n • j 2 = uiM'(u 0 )V n /zi = ~MiVnMi, 


(47) 


which is not in general zero. 


Inner expansion about the interface 

From the 0(1) terms in (13211 , we obtain 

d p (M(U 0 )d pm ) = 0. (48) 

Thus, M(Uo)d p rj 2 is constant in p and since we can identify this expression via (I34p as 
J n ,i, which has to match with n • ji =0. Therefore we can deduce that 

Jn, i = M(Uo)d p r]2 = 0, (49) 

and 7/2 is independent of p. The solution for 7/2 is found in essentially the same way as 
in Section [3l see (fl9ll (l26l) . thus 

V2(s,t) = ^. (50) 
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0(e 3 ) correction 

Noting that r]Q, r]i and 772 are independent of p, the O(e) terms in (1321) yield 

-v n d p U Q = d p M{U 0 )d p r, 3 + ^M(U 0 )d ss K. (51) 

Integrating equation (1511) from — oo to oo, we arrive at 

V n = ^[M(U 0 )d p Tl3]-oo + \ d ssK- (52) 

From (I34p . we can identify the term in the bracket as 

<4,2 = M(U 0 )d p r] 3 . (53) 

At p —Y — oo, we need to match 773 and J n}2 with the solution for 773 and n-J2 in the inner 
layer at F, which in the former case is a function independent of z, and in the latter is 
just zero. Thus, 773 is matched to a constant for p —> — 00 , and J U) 2 is matched to zero, 
thus 

lim M(U 0 )d p r) 3 = lim J n2 = 0. (54) 

p — y —oo p—y— oo 

We next consider the contribution from J n2 as p —» 00 . It is tempting to use (I53|) 
to argue that, since M{Uq) —>■ 0 exponentially fast, J n p also has to tend to zero. Then, 
however, J U)2 cannot be be matched with n • j 2 , as we cannot simply set the latter to 
zero: The bulk equation (|45[) has already got a boundary condition at £, namely ()46l) . 
and setting n • j 2 = 0 would impose too many conditions there. We therefore drop the 
idea that J n p —> 0 as p —> 00 and match the normal fluxes, 

lim J n 2 = n • j 2 b , (55) 

p—yoo s 

Keeping in mind that non-trivial solutions for p\ will arise from (145P , (1461) and V n //i = 
0 at dfl, we expect that J n ^ 2 will not, in general be zero because of (147[) and (l55l) . 
Substituting (153[) and ()47[) into the left and right hand sides of (1551) . respectively, we 
obtain 

lim M(Uo)d p ri 3 = -/riV ny ui| c , (56) 

p —^00 z 

so that now the boundary terms in (152[) have been determined in terms of // 1 . Now, 
however, we have to accept that in general there will be exponential growth in 7/3 as 
p —> 00 : if the left hand side of (I56p is nonzero, and M(Uq) —>• 0 exponentially fast as 
p —>• 00 , then 773 has to grow exponentially. In fact, if we solve (1531) for 773 , and eliminate 
Jn,2 via (|55D and (147p . we obtain the solution 

m = -—- (e p + 2p) + 773 , (57) 
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where 7/3 is an integration constant. The term proportional e 2p is the exponentially 
growing term and it does not appear to be matchable to the outer solution. We will 
resolve this issue in a separate section, by introducing another inner layer, and for now 
continue with analysing the sharp interface model, which in summary is given by 

V • OiV/Ui) = 0, in fi+, 

2 * 

!M = gK, on C, 

//1 — 0, on 5f7 e xt; 

2 1 

Vn = 2 d ss K + -/XlVn/Ul On £ 

4.4 Additional inner layer 

The exponential growth of 773 at p —> 00 is a direct consequence of the exponential decay 
of M(Uq) to 0 as Uq approaches —1 exponentially fast. Notice, however, that the inner 
solution including the correction terms does not decay to —1, because U\(p —» 00 ) > 0, 
so that 

M(U 0 + sU\ + ••■) = M(Uq) + £M'{Uq)U\ + ■■■ 

approaches a non-zero 0(e) value as p —> 00 . We need to ensure that the correction 
eM'(U q)U\ to M(Uq) enters into the calculation of the chemical potential as soon as p 
is in the range where M(Uq) and eM' (Uq)U\ have the same order of magnitude. This 
happens when Uq + 1 = O(e), i.e. when p ~ —(1/2) lne. We therefore introduce another 
layer via 


(58a) 

(58b) 

(58c) 

(58d) 



U(y) 


U(p), fj(y) = Tj(p ), J(y) = J (p). 


Notice the similarity with the change of variables at T. Indeed, the solution in the 
new layer will have exponential terms in the expansion at y —> —00 that need to be 
matched with the expansion at the interface p —> 00 . In terms of the new variables, the 
Cahn-Hilliard equation becomes 


(59) 


(60) 

We expand 

U = —1 + eUi + e^U -2 + • • • , f) = ef/i + s 2 772 + ■ ■ ■ , 

J = Jo + eJ 1 + e 2 J2 + ■ ■ • , 


e 2 d T U - £v n d y U = V • 


(M([>)Vt)) , 




£K 


1 + en (y — | In e) 


-d. 


dyU 


d,U 


1 + £K (y — \ In e) l 1 + £k (y — | In e) 


+ f\U). 
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where we have tacitly anticipated that 7)0 = 0, J_i = 0. Inserting these gives 


V 


= d v 


M'i-lfrdyfk 


+ £dy 


M’ (-l)Uid y f/2 +0(e 2 ). (61) 


The normal flux j n = n • J is given by 

Jn = = [m'(-1)17i + e ((M"(-l)/2) U 2 + M'(-1)U 2 ) + 0(e 2 ) 

x [ edyfji + £ 2 dyf)2 + 0(e 3 )] • 


(62) 


Comparison with the ansatz for the expansion of J immediately implies J n ^ = 0. 


Leading order problem 

To leading order, we have 


-d v 


M'i-yutdyfn 


= 0 , 


— dyyUl + f"(~l)Ul = fjl- 


(63) 


Integrating the first of these once, we obtain that the expression in square bracket has 
to be a constant in y. From (l62l) . we see this is the term J n q in the normal flux, which 
has to match to J n q and n ■ ji in the interface layer and the outer problem, respectively. 
Thus J n \ = 0. Therefore, the contribution 771 is also a constant that needs to match to 
the same value k/6 towards the outer and the interface layer, i.e. for y —> ± 00 , so that 
we have 

m = Ui = c\e~ 2y + c 2 e 2y + (64) 

3 6 

Matching this to the constant outer u\ = k/6 , obtained from (1401) and (l43l) . forces c 2 = 0. 
We next expand Uq at p —> 00 , 

U 0 = -1 + 2e~ 2p + 0(e -4p ). (65) 


The second term accrues a factor of e upon passing to y-variables, and thus has to match 
with the exponential term in eUi, giving c\ = 2 and 


U\ = 2e~ 2y + -k. 

6 


First correction problem 

To next order, we obtain 


dyyU 2 


-dv 


M' {—l)Uid y f}2 


Kd y Ui + f"{-\)U2 + f'"{-\)lh 


J, 


n, 2 


0 , 

m-, 

M'(—l)Uid y fj2- 


( 66 ) 


(67a) 

(67b) 

(67c) 
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£ 

0.01 

0.005 

0.003 

0.002 

0.001 

Eq CSD 

Eq (ED 

CN 

II 

s 

-133.2 

-133.8 

-136.0 

-136.3 

-137.0 

-137.4 

-128 


Table 1: Relaxation rates obtained from the linearised phase field model (17311 are shown 
for different values of £ in the first five columns, and compared to the eigenvalues ob¬ 
tained for linearised sharp interface models for pure surface diffusion (1711) and the porous 
medium type model dza in the next-to-last and the last column, respectively, with 
m = 2/3. 


From (I67al) and (I67cl) . and matching the flux contribution J U) 2 to the outer n • J 2 , we 
obtain 


M'{—l)U\d y fi2 = - /JiV n/ ui| c , 

which in turn has the solution 


m = 


/eM'(-l) 11 V12 


e 2y + 1) + 


36' 


( 68 ) 


(69) 


The integration constant has been fixed by matching r )2 for V —00 with the interface 
solution r/ 2 , see tfUD. We now need to check if the exponential term in (1691) matches with 
the exponential term in (1571) . Expanding at y —> —00 is trivial, and then substituting in 
y = p + In e/2 gives 


m 


8M'(-1) MlVnWl ^ 


2 p + — 


36 


(70) 


Thus, e 2 f )2 contains a term proportional to e 3 e 2?/ term that is identical to the e 3 e 2?/ term 
that appears in e 3 %, see (1571) . Thus, we have resolved the issue with the exponentially 
growing term (for p —> 00 ) in the correction to the chemical potential in the interface 
layer expansion. 


4.5 Linear stability analysis 

Besides the usual surface diffusion term, equation (158|1 contains an additional normal 
flux term which is nonlocal. In cases where there are multiple regions of u close to 1, the 
nonlocal term couples the interfaces of these regions with each other and drive coarsening 
where the larger regions grow at the expense of smaller ones. This is not expected for 
pure surface diffusion. Even for a single convex domain that is slightly perturbed from 
its radially symmetric state, the effect on the relaxation dynamics is noticeable, as we 
now explore. 

To compare the sharp interface model with the phase field model, we consider the 
relaxation of an azimuthal perturbation to a radially symmetric stationary state with 
curvature k = 1/ro- For azimuthal perturbations proportional to cos md, the pure surface 
diffusion model v n = yjld ss K predicts an exponential decay rate 

a = 11 . 

Tn 
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In contrast, the decay rate in the porous medium model, Equation ([551) . is given by 


A 


2 m 2 (m 2 — 1) 

3 


1 m(m 2 — 1) 
9 ^ 


tanh(mlogr 0 1 ). 


In the diffuse interface model, the perturbation v\ (r, t) cos m0 satisfies 


vit 


mi 


\ d ( 

-ref ( rM( ”“ ) 
e 2 d ( dv\ 
r dr \ dr 


9mi 

dr 



) ~ ^-M(u 0 )mi, 
(“) v i + f"( v o)vi 


(72) 


(73) 


where vo(r) is the radially symmetric stationary state. We solve this system numerically, 
using the Chebyshev spectral collocation method (see Appendix) with At = 10~ 3 and 400 
mesh points until t = 1/e 2 . The decay rate of the eigenfunction is tracked by monitoring 
its maximum. The diffuse interface decay rates are scaled with 1 /e 2 to compare with the 
sharp interface model. The base state that is needed for this calculation is determined 
a priori with the interface, i.e. the zero contour, positioned at ro = 0.5. The initial 
condition for the perturbation, 

wi(0,r) = exp [l/(a 2 - (r 0 - r) 2 )], (74) 

acts approximately as a shift to the leading order shape of the inner layer. The constant 
a is chosen so that the support of ui(0,r) lies in the range r > r*. 

The results are compared in Table [T} They show that the decay rate of the azimuthal 
perturbation to the radially symmetric base state obtained for m = 2 tends to the 
eigenvalue for the linearised sharp interface model with the contribution from nonlinear 
bulk diffusion, rather than to the one for pure surface diffusion. This confirms that 
(1581) describes the leading oder sharp interface evolution for the Cahn-Hilliard model Jl]) 
correctly, and that the sharp interface motion is distinct from the one induced by pure 
surface diffusion. 


5 Modifications 

5.1 Solutions with u > 1 for the mobility M(u) = |1 — u 2 \ 

As pointed out in Section [3j solutions that have a modulus |rt| > 1 and converge to the 
usual stationary Cahn-Hilliard solutions are conceivable for the mobility M(u ) = |1 — u 2 \ 
and are seen to arise in numerical solutions with this mobility for appropriate initial 
conditions. For this case, we can carry out the asymptotic derivations to obtain the 
sharp interface limit and match the inner problem to outer solutions on both sides of 
the interface, accepting thereby that the outer solution for u in is larger than one. 
Otherwise the detailed derivations follow the same pattern as in section l~T~3l and can be 
found in |40j . 
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£ 

0.01 

0.005 

0.002 

0.001 

Eq dZHD 

\n =2 

-144.7 

-146.3 

-147.5 

-147.8 

-148.1 


Table 2: The decay rates of an azimuthal perturbation obtained by the diffuse and sharp 
interface models show good agreement for general initial condition not bounded between 
±1 and mobility M(u ) = 1 — u 2 . The numerical method and discretisation parameters 
are the same as in Table [T] The description of the numerical approach and parameters 
carries over from Table Q] 


The upshot is that the sharp interface model now has contributions from nonlinear 
bulk diffusion on both sides of the interface, in addition to surface diffusion, viz. 


V • (nf'Vnf) = 0, on Q±, 

± 2 

Hi = -k, on C, 

V„/i| = 0, on <9fl, 

2 1 

v n — —dssK T ^(/T ^nP\ T AT ^n/T )j on C- 

This sharp interface model predicts an exponential decay rate of 
2 m 2 (m 2 — 1) 1 m(m 2 — 1) 


A = — 


-(tanh(mlogr 0 x ) + 1) 
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(75a) 

(75b) 

(75c) 

(75d) 


(76) 


for the evolution of the perturbation to the radially symmetric stationary state with wave 
number m. Table [2] shows that equation (17611 is indeed consistent with numerical results 
for the diffuse model. As a cautionary remark, we note that we are dealing here with 
a sign-changing solution of a degenerate fourth order problem, in the sense that 1 — u 
changes sign and the mobility degenerates. The theory for this type of problems is still 
being developed mmmmmm- 


5.2 Degenerate biquadratic mobility 

For the mobilities investigated so far, nonlinear bulk diffusion enters at the same order 
as surface diffusion. If we employ M(u ) = ((1 — u 2 ).|_) 2 , then 

32 = uiM\u 0 )\7 n Hi = 0. (77) 

The contribution of the bulk diffusion flux to the normal velocity of the interface is 
subdominant to surface diffusion and therefore 

1 4 

v n = - sech 4 p dp d ss K = -d ss K. (78) 

^ J —OO ^ 

Table [3] shows that the decay rate obtained from the numerical solution of the diffuse 
interface model for the degenerate biquadratic mobility is indeed consistent with the 
predictions obtained for the sharp interface model (17811 with pure surface diffusion. 
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£ 

0.01 

0.005 

0.001 

Eq m 

A m —2 

-84.6 

-84.7 

-85.2 

-85.3 


Table 3: The decay rates obtained by the diffuse interface model for the mobility M(u ) = 
((1 — u 2 ) + ) 2 and |u| < 1 show good agreement with the surface diffusion model in (jTTI) . 
with 9JT = 4/9, as e —>• 0. The description of the numerical approach and parameters 
carries over from table (fTj). 

6 Conclusions 

In this paper, we have derived the sharp interface limit for a Cahn-Hilliard model in 
two space dimensions with a nonlinear mobility M(u) = (1 — 'it 2 )_|_, and a double-well 
potential with minima at ±1 for the homogeneous part of the free energy. We found that 
in addition to surface diffusion, there is also a contri bution from bulk diffusion to the 
interface motion which enters at the same order. This contribution enters only from one 
side of the interface, whereas for the mobility M(u ) = |1 — u 2 |, solutions have also been 
considered for which bulk diffusion in the sharp interface limit enters from both sides at 
the same order as surface diffusion. 

The situation studied here was focused on the case of convex = {x G 17; u > 0} 
with an 0(1) curvature for the interace u = 0, though the asymptotic analysis also 
remains valid if fl_|_ is the union of well-separated convex domains. The dynamics for 
concentric circles of different phases has also been looked into m ■ For the case where the 
interface has turning points, the derivation needs to be revisited, since the the location 
of the free boundary T, given by p = er in inner coordinates about the interface, depends 
on the curvature so that |cr| —> oc if k tends to zero. Moreover, as the curvature changes 
sign, T changes the side of the interface. On a different plane, it would also be interesting 
to investigate the coarsening behaviour m for the sharp interface model (15811 . For 
ensembles of two or more disconnected spheres, pure surface diffusion does not give rise 
to coarsening, but coarsening is expected for the mixed surface/bulk diffusion flux in 

m- 

While the Cahn-Hilliard equation © plays a role in some biological models, see for 
example m, and may have significance in modelling spinodal decomposition in porous 
media, possibly with different combinations of mobilities, e.g. M(u) = |1 — u 2 \ + a( 1 — 
u 2 ) 2 , see [40] , the main motiviation for our investigation stems from the role degenerate 
Cahn-Hilliard models play as a basis for numerical simulations for surface diffusion with 
interface motion driven by (J2]). The upshot for the specific combination of mobility and 
double well potential used in dTJ) is not useful for this purpose, since a contribution from 
bulk diffusion enters at the same order. For mobilities with higher degeneracy, such as 
M(u) = ((1 — u 2 )+) 2 , this undesired effect is of higher order and can be made arbitrarily 
small, at least in principle, by reducing e. Nevertheless, for finite e. it is still present and 
a cumulative effect may arise for example through a small but persistent coarsening of 
phase-separated domains. 

A range of alternatives can be found in the literature, in particular using the combi- 
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nation of M = (1 — w 2 )+ or M = |1 — u 2 \ with the logarithmic or with the double obstacle 
potential |18| . These combinations force the order parameter u to be equal to or much 
closer to ±1 away from the interface, thus shutting out the bulk diffusion more effec¬ 
tively. Numerical methods have been developed for these combinations and investigated 
in the literature, see for example [3 00® METIS. Other approaches that have been 
suggested include a dependence of the mobility on the gradients of the order parameter 
[42] . tensorial mobilities 133] . or singular expressions for the chemical potential [50] . 

As a final remark, we note that many analytical questions remain open. For example, 
the existence of solutions that preserve the property that |u| > 1 in some parts of U has 
not been shown. Also, the approximation or (0) by a free boundary problem Q should 
be investigated systematically using b = min (1 — |u|) > 0 as a small parameter, in 
the spirit of what was done, for example, in [35] for the precursor model of a spreading 
droplet. The conditions at the free boundary T could then be recovered from matching 
to an inner solution. If fo —>• 0 in finite time, the effect of the “precursor” regularisation 
is lost and either the regularising effect implicit in the numerical discretisation or any 
explicit regularisation that is used (e.g., the one suggested in [23]) have to be taken into 
account. It would be interesting to see for which regularisations the conditions in (l3cl) 
are recovered. We note, however, that the evolution of the leading order sharp interface 
model in is usually insensitive to the conditions imposed at T. 

7 Appendix: Numerical Methods 

We numerically solved the radially symmetric counterpart to (0) in polar coordinates 
without an explicit regularsisation (such as the on used in [23]) via a Chebyshev spectral 
collocation method in space and semi-implicit time-stepping, using a linearised convex 
splitting scheme to treat /. For details on spectral methods in general, we refer the 
reader to the references [551 [56] . We also split the mobility as M (u) = (M(u) — 9) + 9, to 
evaluate ( M(u ) — 6) at the previous time step whilst solving the remaining 9 portion at 
the next time step, which improved the stability. We choose 9 = O.Ole in our simulations. 
Varying 9 confirmed that the results did not sensitively depend on its value provided it 
was O(s). 

As the Chebyshev-Lobatto points are scarcest in the middle of the domain, we resolve 
the interior layer by introducing a non-linear map x £ [—1,1] i—>■ r £ [0,1], as suggested in 
[14] . r = (1/2) + arctan (<5tan7rx/2) /7T, where 0 < 5 < 1 is a parameter that determines 
the degree of stretching of the interior domain, with a smaller value of 5 corresponding 
to greater degree of localisation of mesh points about the centre of the domain. In this 
paper, we general set 5 = lOe. This choice of 5 is guided by numerical experiments, which 
show that further increase in the number of mesh points does not alter the stationary 
solution. Moreover, since r = 0 is a regular singular point, we additionally map the 
domain r € [0,1] linearly onto a truncated domain [10~ 10 ,1]. Again, we verified that 
varying the truncation parameter did not affect the numerical results. Unless otherwise 
stated, the numerical simulations reported in the paper are done with 400 collocation 
points and timestep At = 10~ 3 . 
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The linearised phase-field models were solved using the same method, with a base 
state that was obtained from a preceding run and then “frozen” in time, i.e. not co-evolved 
with the perturbation. 
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